Exercício computacional I

Dataframes podem ser interpretados como tabelas de dados ou dados tabulados e, talvez, sejam uma das estruturas de dados linha-coluna mais importantes na linguagem R no momento em que pensamos em aplicá-la para ciência de dados. Por isso, esse exercício tem um propósito básico: criar e explorar comandos básicos relacionados com dataframes no R.

  1. Crie o dataframe mostrado na figura acima e armazene no objeto df
df <- data.frame(
  id = c (1:5),
  Empresa = c("A","B","C","D","E"),
  Indices = c(500.3,530.2,630.5,400.20,940.20),
  Datas = as.Date(c("2020-03-05", "2020-04-21", "2020-12-10", "2020-10-15","2020-09-20")),
  stringsAsFactors = FALSE)
df
##   id Empresa Indices      Datas
## 1  1       A   500.3 2020-03-05
## 2  2       B   530.2 2020-04-21
## 3  3       C   630.5 2020-12-10
## 4  4       D   400.2 2020-10-15
## 5  5       E   940.2 2020-09-20
  1. Utilize a função str() e interprete os resultados sobre cada tipo de dado contido no dataframe
str(df)
## 'data.frame':    5 obs. of  4 variables:
##  $ id     : int  1 2 3 4 5
##  $ Empresa: chr  "A" "B" "C" "D" ...
##  $ Indices: num  500 530 630 400 940
##  $ Datas  : Date, format: "2020-03-05" "2020-04-21" ...

Neste dataframe temos 5 observações e 4 variáveis. As variáveis encontradas são id (int), Empresa (chr), Indices (num), Datas (Date)

  1. Faça a extração apenas das colunas de empresas e índices
df2 <- data.frame(df$Empresa, df$Indices)
df2
##   df.Empresa df.Indices
## 1          A      500.3
## 2          B      530.2
## 3          C      630.5
## 4          D      400.2
## 5          E      940.2
  1. Crie um array com os elementos relacionados com: a primeira (1) e terceira (3) linhas e a segunda (2) e quarta (4) colunas.
array_elementos <- df[c(1,3),c(2,4)]
array_elementos
##   Empresa      Datas
## 1       A 2020-03-05
## 3       C 2020-12-10
  1. Adicione uma nova coluna ao dataframe com os setores empresariais “IT”, “adm”, “executivo”, “RH”, “O&M” e armazene em novo dataframe chamado df3
df3 <- df
df3$setores <- c("IT", "adm", "executivo", "RH", "O&M")
df3
##   id Empresa Indices      Datas   setores
## 1  1       A   500.3 2020-03-05        IT
## 2  2       B   530.2 2020-04-21       adm
## 3  3       C   630.5 2020-12-10 executivo
## 4  4       D   400.2 2020-10-15        RH
## 5  5       E   940.2 2020-09-20       O&M
  1. Combine o dataframe do item 1), dado por df, com o novo dataframe mostrado abaixo e armazene o resultado, também como dataframe, no objeto dfn. Estude as funções rbind e cbind para isso.

dfn <- data.frame(
  id = c(6:10),
  Empresa = c("F", "F", "G", "K", "L"),
  Indices = c(1200.3, 230.4, 100.5, 905.4, 1100.5),
  Datas = c("2020-09-10", "2020-07-08", "2020-10-15", "2020-06-07", "2020-02-22"),
  stringsAsFactors = FALSE
)

dfn
##   id Empresa Indices      Datas
## 1  6       F  1200.3 2020-09-10
## 2  7       F   230.4 2020-07-08
## 3  8       G   100.5 2020-10-15
## 4  9       K   905.4 2020-06-07
## 5 10       L  1100.5 2020-02-22

dfn <- rbind(df, dfn)
dfn
##    id Empresa Indices      Datas
## 1   1       A   500.3 2020-03-05
## 2   2       B   530.2 2020-04-21
## 3   3       C   630.5 2020-12-10
## 4   4       D   400.2 2020-10-15
## 5   5       E   940.2 2020-09-20
## 6   6       F  1200.3 2020-09-10
## 7   7       F   230.4 2020-07-08
## 8   8       G   100.5 2020-10-15
## 9   9       K   905.4 2020-06-07
## 10 10       L  1100.5 2020-02-22

Exercício Computacional II

Conceito explorado: Geração de Números Aleatórios

Gere uma amostra com 1000 observações que segue a distribuição de probabilidade Gaussiana com média μ = 10 e desvio padrão σ = 5, Armazene os números aleatórios gerados no objeto r.

?norm

rnorm é a função que simula variáveis normais. A forma de usar é:

x <- rnorm(n, mean, sd)

em que n é o tamanho da amostra (x será um vetor caso n > 1), e mean e sd são parâmetros (opcionais) dando a média e o desvio-padrão da normal. Se mean ou sd forem omitidos, serão usados os valores, respectivamente, de 0 e 1.

n      = 1000
media  = 10
desvio = 5
r <- rnorm(n, media, desvio)
head(r)
## [1] 16.779313 18.634569  4.384529 13.213542  9.313709  4.926157
  1. Qual é o tipo de objeto r? Quais instruções você utilizou para verificar essa informação?
str(r)
##  num [1:1000] 16.78 18.63 4.38 13.21 9.31 ...
summary(r)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -8.122   6.469  10.202   9.944  13.629  29.436
  1. Obtenha o histograma relacionado com o vetor r.
hist(r, probability = TRUE, ylim = c(0,0.1), breaks = 20, main = 'Estimação do Histograma')

  1. Plote, sobre o histograma, a curva de densidade normal informando os valores de média e desvio padrão. Dica: no R, as funções curve e dnorm são úteis para solucionar esse ponto.
?dnorm
# Esta função está dando erro ao ser gerada no compiler do pdf, verificar.
#curve(dnorm(x, media, desvio), add = TRUE, col = 'red', lwd = 1)
# dnorm -> d é densidade
  1. Utilize o pacote ggplot2 da linguagem R para obter o mesmo resultado dos itens anteriores.
ggplot(data.frame(r), aes(x = r)) + 
  geom_histogram(aes(y =..density..),
                 bins   = 20, 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm, args = list(mean = media, sd = desvio))

?stat_function

Exercício Computacional III

Conceito explorado: Estatística Descritiva e Análise Exploratória de Dados

A Figura abaixo mostra um instrumento de teste (Field Fox Keysight) que pode ser usado em labo- ratório ou em campo para medições de sinais de radiofrequência como os presentes em sistemas de comunicações sem fio. Isso significa que podemos usar esse equipamento para análise de redes sem fio, cobertura de operadoras de telecomunicações, além de testes com dispositivos de RF e outros equipamentos de telecomunicações. Nesse contexto, utilizamos esse equipamento para a realização de medições de intensidade de sinal no campus do Inatel a fim de levantar a cobertura de uma rede sem fio experimental, configurada para transmitir sinais na faixa de frequência de ondas milimétricas. O estudo de cobertura e propagação nessa faixa de frequência é um aspecto de pesquisa relevante para sistemas de comunicações da quinta geração de redes móveis. Nesse exercício, temos o objetivo de fazer a análise exploratória de dois conjuntos de dados, dataset_1 e dataset_2 exportados pelo instrumento de teste.

  1. Faça a importação do arquivos dataset_1 e dataset_2 exportados pelo equipamento para o ambiente do RStudio.
dataset_1 <- read.delim("Datasets/dataset_1.csv",stringsAsFactors = FALSE)
dataset_2 <- read.delim("Datasets/dataset_2.csv",stringsAsFactors = FALSE)
  1. Análise o resultado da importação, como as estruturas e tipos de variáveis. Quais são as principais informações contidas no arquivo?
head(dataset_1, 21)
##                                          X..FILETYPE.CSV
## 1                                        ! VERSION 1.0,1
## 2             ! TIMESTAMP Sunday, 08 April 2018 17:54:25
## 3                        ! TIMEZONE (UTC-03:00) Brasilia
## 4                           ! NAME Keysight Technologies
## 5                                         ! MODEL N9952A
## 6                                    ! SERIAL US55240250
## 7                             ! FIRMWARE_VERSION A.09.12
## 8                                          ! CORRECTION 
## 9                                       ! Application SA
## 10               ! Trace TIMESTAMP: 2018-04-08 17:54:25Z
## 11                                   ! Trace GPS Info...
## 12                         ! GPS Latitude: 22 15.40727 S
## 13                        ! GPS Longitude: 45 41.74932 W
## 14                      ! GPS Seconds Since Last Read: 0
## 15                                 ! CHECKSUM 0371155866
## 16 ! DATA Freq,SA Clear-Write,SA Blank,SA Blank,SA Blank
## 17                                        ! FREQ UNIT Hz
## 18                                       ! DATA UNIT dBm
## 19                                                 BEGIN
## 20                   28000002997,-88.5112531089517,0,0,0
## 21                   28000002997,-88.4285038898773,0,0,0
head(dataset_2, 21)
##                                          X..FILETYPE.CSV
## 1                                        ! VERSION 1.0,1
## 2           ! TIMESTAMP Sunday, 01 January 2006 13:47:50
## 3      ! TIMEZONE (UTC-08:00) Pacific Time (US & Canada)
## 4                           ! NAME Keysight Technologies
## 5                                         ! MODEL N9952A
## 6                                    ! SERIAL US55240250
## 7                             ! FIRMWARE_VERSION A.09.12
## 8                                          ! CORRECTION 
## 9                                       ! Application SA
## 10               ! Trace TIMESTAMP: 2006-01-01 13:47:50Z
## 11                                   ! Trace GPS Info...
## 12                        ! GPS Latitude: 22 15.382499 S
## 13                       ! GPS Longitude: 45 41.758792 W
## 14                      ! GPS Seconds Since Last Read: 0
## 15                                 ! CHECKSUM 0621167751
## 16 ! DATA Freq,SA Clear-Write,SA Blank,SA Blank,SA Blank
## 17                                        ! FREQ UNIT Hz
## 18                                       ! DATA UNIT dBm
## 19                                                 BEGIN
## 20                    28000004352,-106.67464776889,0,0,0
## 21                     28000004352,-106.6737572421,0,0,0

Mediante a observação dos dados carregados é possível dizer que os dados a serem utilizados estão conditos a partir da linha 21.

Criando função para extrair informações das medidas.

obter_medidas <- function(data) {
   data[21:nrow(data)-1,]
}

Obtendo informações das medidades do dataset_1

medidas_1 <- obter_medidas(dataset_1)
head(medidas_1)
## [1] "28000002997,-88.5112531089517,0,0,0" "28000002997,-88.4285038898773,0,0,0"
## [3] "28000002997,-88.6968101484764,0,0,0" "28000002997,-88.1224479155788,0,0,0"
## [5] "28000002997,-88.090961444615,0,0,0"  "28000002997,-87.6165186878635,0,0,0"

Obtendo informações das medidades do dataset_1

medidas_2 <- obter_medidas(dataset_2)
head(medidas_2)
## [1] "28000004352,-106.67464776889,0,0,0"  "28000004352,-106.6737572421,0,0,0"  
## [3] "28000004352,-106.672865298382,0,0,0" "28000004352,-106.671971938544,0,0,0"
## [5] "28000004352,-106.671077163395,0,0,0" "28000004352,-106.670180973743,0,0,0"

Cria função para extrair amostrar de coleta

# A partir de uma amostra descobre onde se encontra a primeira "," 
amostra <- medidas_1[21]
local <- str_locate(amostra,",")

create_amostras_df = function(medidas, local) {
  data <- substr(medidas, start = local[1]+1, stop = local[1]+10) %>% 
    as.numeric() %>% data.frame()
  colnames(data) <- "Prx"
  data
}
df1_table <- create_amostras_df(medidas_1, local)
head(df1_table)
##         Prx
## 1 -88.51125
## 2 -88.42850
## 3 -88.69681
## 4 -88.12245
## 5 -88.09096
## 6 -87.61652
df2_table <- create_amostras_df(medidas_2, local)
head(df2_table)
##         Prx
## 1 -106.6746
## 2 -106.6737
## 3 -106.6729
## 4 -106.6720
## 5 -106.6711
## 6 -106.6702
  1. Obtenha o histograma dos valores de potência de recepção coletados pelo equipamento em cada conjunto de dados.
hist(df1_table$Prx, probability = TRUE, breaks = 100, main = 'Estimação do Histograma df_table1')

# Apresentando duas formas de gerar o histograma.
ggplot(df2_table, aes(x = Prx)) + geom_histogram(aes(y =..density..),
                                                 bins   = 50, 
                                                 colour = "black", 
                                                 fill = "white")

hist(df2_table$Prx, probability = TRUE, breaks = 100, main = 'Estimação do Histograma')

  1. Em qual localidade específica foram realizadas as medições de cada conjunto de dados?
get_coords <- function(data) {
  gps_lat_ponto_1  = data[12,1]
  # gps_lat_ponto_1
  gps_long_ponto_1 = data[13,1]
  # gps_long_ponto_1
  
  # Conversão de grau, minuto fracionado
  grau_lat_ponto_1        = substr(gps_lat_ponto_1, start = 17, stop = 18)
  minuto_lat_ponto_1      = substr(gps_lat_ponto_1, start = 19, stop = 26)
  grau_long_ponto_1       = substr(gps_long_ponto_1, start = 18, stop = 19)
  minuto_long_ponto_1     = substr(gps_long_ponto_1, start = 20, stop = 27)
  
  # Composição dos valores numéricos de latitude e longitude
  data.frame(
    lat = (-1)*(as.numeric(grau_lat_ponto_1) + as.numeric(minuto_lat_ponto_1)/60),
    lng = (-1)*(as.numeric(grau_long_ponto_1) + as.numeric(minuto_long_ponto_1)/60)
  )
}

coord <- get_coords(dataset_1)
coord
##         lat       lng
## 1 -22.25679 -45.69582

Plotando as coordenadas

leaflet() %>%
  addTiles() %>%
  addMarkers(lng = coord$lng, lat = coord$lat,
             popup = 'ponto a',
             clusterOptions = markerClusterOptions())
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Exercício Computacional IV

Conceito explorado: Estatística Descritiva e Análise Exploratória de Dados

Considere o mesmo contexto do exercício anterior e um conjunto maior de arquivos .csv, que são exportados pelo instrumento de medida e armazenados em um diretório.

  1. Com os arquivos .csv armazenados no diretório, elabore uma rotina em linguagem R para fazer a leitura de todos os arquivos de forma otimizada.

Lendo dados dos arquivos e adicionando o prefixo do path.

temp_files <- paste0("Datasets/", list.files(path = "Datasets", pattern = ".csv"))
temp_files
##  [1] "Datasets/dataset_1.csv"  "Datasets/dataset_10.csv"
##  [3] "Datasets/dataset_2.csv"  "Datasets/dataset_3.csv" 
##  [5] "Datasets/dataset_4.csv"  "Datasets/dataset_5.csv" 
##  [7] "Datasets/dataset_6.csv"  "Datasets/dataset_7.csv" 
##  [9] "Datasets/dataset_8.csv"  "Datasets/dataset_9.csv"

Unificando os dados do dataset em um dataset único

files <- lapply(temp_files, read.delim)
  1. Capture os dados de geolocalização (latitude, longitude) de todos os arquivos, faça os processamentos e transformações necessárias, visando o armazenamento em um dataframe.
get_coords <- function(data) {
  gps_lat_ponto_1  = data[12,1]
  # gps_lat_ponto_1
  gps_long_ponto_1 = data[13,1]
  # gps_long_ponto_1
  
  # Conversão de grau, minuto fracionado
  grau_lat_ponto_1        = substr(gps_lat_ponto_1, start = 17, stop = 18)
  minuto_lat_ponto_1      = substr(gps_lat_ponto_1, start = 19, stop = 26)
  grau_long_ponto_1       = substr(gps_long_ponto_1, start = 18, stop = 19)
  minuto_long_ponto_1     = substr(gps_long_ponto_1, start = 20, stop = 27)
  
  # Composição dos valores numéricos de latitude e longitude
  data.frame(
    lat = (-1)*(as.numeric(grau_lat_ponto_1) + as.numeric(minuto_lat_ponto_1)/60),
    lng = (-1)*(as.numeric(grau_long_ponto_1) + as.numeric(minuto_long_ponto_1)/60)
  )
}

# Inicializando os vetores
lat  = vector(mode="numeric", length = length(files))
lng = vector(mode="numeric", length = length(files))

for(i in 1:length(files)) {
  coord <-  get_coords(files[[i]]);
  lat[i] <- coord$lat
  lng[i] <- coord$lng
}

# Criando o dataframe de coordenadas.

coords <- data.frame(lat, lng) %>% mutate(Pontos = c("Ponto 1", "Ponto 2", "Ponto 3", "Ponto 4", "Ponto 5", "Ponto 6", "Ponto 7", "Ponto 8", "Ponto 9", "Ponto 10"))
  1. Apresente no mapa os dados de geolocalização obtidos no item anterior.
leaflet() %>% 
  addTiles() %>% 
  addMarkers(lng = coords$lng, lat = coords$lat, label = coord$Pontos,
             popup = paste0(coords$Pontos, "<br/>Lat: ", coords$lat, "<br/>Lng: ", coords$lng))

Exercício Computacional V

Conceito explorado: MSE - Mean Square Error

Considere o seguinte modelo de geração de dados mostrado abaixo: \[ y = h(x) +ε \]

(1.3)

Nesse modelo, h(x) = 3x+30 consiste na função hipótese verdadeira, muitas vezes desconhecida na prática de ML, e ε é um termo que expressa a incerteza entre os valores da função hipótese verdadeira e a variável de saída ou resposta y. Estatisticamente, ε é interpretado como um ruído, que nesse exercício segue a distribuição de probabilidade Gaussiana com média μ = 0 e desvio padrão σ = 15. A notação em negrito usada ocorre em função de (1.3) ser um modelo vetorial de dados. A variável explanatória x usada será um vetor de valores inteiros de zero (1) a cem (100). Logo, o modelo de dados é formado pelos vetores y, h(x) e ε, sendo cada um com dimensões (número de linhas e colunas) de 100×1. Considere que um grupo de cientistas de dados já realizaram o trabalho de modelagem e encontraram uma função hipótese candidata dada por: \[ hˆ(x) = 2.8x+32 \]

(1.4)

  1. Construa esse modelo de geração de dados. Para que seja possível a reprodução de resultados em função do vetor aleatório ε utilize a semente (seed) 123 em seu código.
# Quantidade de amostras
n <- 100

#Gerando as amostras
x <- seq(1, 100, length.out = n)

# Função hipótese verdadeira
h_x <- 3*x + 30

set.seed(123)

#Ruido
?rnorm
epsilon <- rnorm(n, 0, 15)

# Variável de saída
y <- h_x + epsilon
  1. Faça um gráfico de dispersão da variável explanatória x com saída conhecida y.
plot(x, y, col=1, pch=1, main = 'Gráfico de Dispersão', 
     col.main = 'black', 
     xlab = 'Variável Explanatória',
     ylab = 'Variável de Saída')
grid()

  1. Obtenha o histograma relacionado com a variável de saída y.
hist(y, main = "Histograma da variável de saída (y)")

  1. A equação do MSE, mostrada abaixo, é uma métrica de desempenho relacionada com qual tipo de tarefa de aprendizagem de máquina? MSE (1.5)

Statistical learning

  1. Faça a estimação do erro quadrático médio do modelo proposto pelos cientistas.
# Função hipótese Estimada
h_x_estimado <- 2.8*x + 32
MSE = (1/n)*sum(((y - h_x_estimado)^2))
MSE
## [1] 320.9044
  1. Faça uma análise: o modelo proposto é plausível para explicar os dados? De quais fatores esse desempenho depende?

** Não ficou muito claro para mim, revisar.

Exercício Computacional VI

Conceito explorado: MSE - Mean Square Error

Considere o mesmo modelo de geração de dados do exercício anterior. O objetivo aqui é constatar o impacto do desvio padrão σ sobre a performance do modelo proposto pelos cientistas de dados. Para isso, utilize a instrução abaixo, em linguagem R, para a geração de um vetor com diversos valores de desvio padrão para a incerteza Gaussiana retratada pelo termo ε. Vetor com valores de desvio padrão entre 0 e 20 std_vector = seq(1, 20, length.out = 100) Para que seja possível explorar e visualizar o impacto de σ - realize, pelo menos, 1000 iterações do algoritmo. Especificamente, para cada valor de desvio padrão avaliado, armazene e faça o cálculo da média aritmética sobre 1000 valores de performance expressos pelo MSE. Um dica é utilizar estruturas em loop (for) para a implementação das iterações.

  1. Construa o modelo de geração de dados incluindo as iterações para cada valor de σ
n = 100
num_iter = 1000

#Construindo o vetor com os valores ordenados em ordem descrescente.
std_vector = seq(1, 20, length.out = 100) %>% sort(decreasing = TRUE)

x = seq(1, 100, length.out = n)

#Inicializando os vetores
MSE_vector = rep(0,length(std_vector))
MSE        = rep(0,num_iter)

std_vector
##   [1] 20.000000 19.808081 19.616162 19.424242 19.232323 19.040404 18.848485
##   [8] 18.656566 18.464646 18.272727 18.080808 17.888889 17.696970 17.505051
##  [15] 17.313131 17.121212 16.929293 16.737374 16.545455 16.353535 16.161616
##  [22] 15.969697 15.777778 15.585859 15.393939 15.202020 15.010101 14.818182
##  [29] 14.626263 14.434343 14.242424 14.050505 13.858586 13.666667 13.474747
##  [36] 13.282828 13.090909 12.898990 12.707071 12.515152 12.323232 12.131313
##  [43] 11.939394 11.747475 11.555556 11.363636 11.171717 10.979798 10.787879
##  [50] 10.595960 10.404040 10.212121 10.020202  9.828283  9.636364  9.444444
##  [57]  9.252525  9.060606  8.868687  8.676768  8.484848  8.292929  8.101010
##  [64]  7.909091  7.717172  7.525253  7.333333  7.141414  6.949495  6.757576
##  [71]  6.565657  6.373737  6.181818  5.989899  5.797980  5.606061  5.414141
##  [78]  5.222222  5.030303  4.838384  4.646465  4.454545  4.262626  4.070707
##  [85]  3.878788  3.686869  3.494949  3.303030  3.111111  2.919192  2.727273
##  [92]  2.535354  2.343434  2.151515  1.959596  1.767677  1.575758  1.383838
##  [99]  1.191919  1.000000
# Função hipótese verdadeira
h_x <- 3*x + 30
h_x
##   [1]  33  36  39  42  45  48  51  54  57  60  63  66  69  72  75  78  81  84
##  [19]  87  90  93  96  99 102 105 108 111 114 117 120 123 126 129 132 135 138
##  [37] 141 144 147 150 153 156 159 162 165 168 171 174 177 180 183 186 189 192
##  [55] 195 198 201 204 207 210 213 216 219 222 225 228 231 234 237 240 243 246
##  [73] 249 252 255 258 261 264 267 270 273 276 279 282 285 288 291 294 297 300
##  [91] 303 306 309 312 315 318 321 324 327 330
for(k in 1:length(std_vector)) {
  for(i in 1:num_iter) {
     # Ruído
    epsilon <- rnorm(n, 0, std_vector[k])
    
    # Variável de saída
    y <- h_x + epsilon
    
    # Função hipótese Estimada
    h_x_estimado <- 2.8*x + 32
    
    # Erro quadrático médio para uma iteração
    MSE[i] = (1/n)*sum(((y - h_x_estimado)^2))
  }
  
  
  # Erro quadrático médio para várias iterações  
  MSE_vector[k] = mean(MSE)
}
  1. Faça um gráfico que mostra o impacto de σ, colocado sobre o eixo x, sobre o desempenho indicado pelo MSE, apresentado no eixo y.
plot(std_vector,MSE_vector,col=1, pch=1, main = 'MSE - Mean Squared Error', 
     col.main = 'black', 
     xlab = 'Desvio padrão dos Erros',
     ylab = 'MSE')
grid()

  1. Faça uma análise: o impacto com o aumento ou redução de σ é significativo para o modelo? Qual a justificativa?

O impacto é sim signficativo, tendo em vista que o MSE tende a crescer com níveis mais altos de ruído.

## ── Attaching packages ────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble  3.0.1     ✓ purrr   0.3.4
## ✓ tidyr   1.0.2     ✓ forcats 0.5.0
## ✓ readr   1.3.1
## ── Conflicts ───────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x purrr::map()    masks maps::map()

Exercício Computacional VII

Conceito explorado: MSE - Mean Square Error

Considere o mesmo modelo de geração de dados do exercício anterior. Agora, nosso objetivo é constatar o impacto do número de amostras n sobre a performance do modelo proposto pelos cientistas de dados. De forma similar ao caso anterior, utilize a instrução abaixo, em linguagem R, para a geração de um vetor com diversos valores de desvio padrão para a incerteza Gaussiana retratada pelo termo ε. Vetor com valores númericos da amostra. n_vector = seq(10, 100, 5)

Para que seja possível explorar e visualizar o impacto de n - realize, pelo menos, 1000 iterações do algoritmo. Especificamente, para cada valor do número de amostras avaliado, armazene e faça o cálculo da média aritmética sobre 1000 valores de performance, expressos pelo MSE. Um dica é aproveitar as estruturas em loop (for) do exercício anterior para a implementação das iterações.

  1. Construa o modelo de geração de dados incluindo as iterações para cada valor de n
num_iter <- 1000
n_vector <- seq(10, 100, 5)

# Inicialização dos vetores de MSE
MSE_vector = rep(0,length(n_vector))
MSE        = rep(0,num_iter)

n_vector
##  [1]  10  15  20  25  30  35  40  45  50  55  60  65  70  75  80  85  90  95 100
for (k in 1:length(n_vector)){
  
  # Crie um vetor de índices x variando de 0 a 100. Utilize a função seq - sequence
  # x = seq(from = 0, to = 100, by = 1)
  n = n_vector[k]
  x = seq(1, 100, length.out = n)
  
  # Função hipótese verdadeira
  h_x <- 3*x + 30
  
  for (i in 1:num_iter){
    
    # Ruído
    std = 0.5
    epsilon <- rnorm(n, 0, std)
    
    # Variável de saída
    y <- h_x + epsilon
    
    # Função hipótese Estimada
    h_x_estimado <- 2.8*x + 32
    
    # Erro quadrático médio para uma iteração
    MSE[i] = (1/n)*sum(((y - h_x_estimado)^2))
  }
  
  
  # Erro quadrático médio para várias iterações  
  MSE_vector[k] = mean(MSE)
}
  1. Faça um gráfico que mostra o impacto de n, colocado sobre o eixo x, sobre o desempenho indicado pelo MSE, apresentado no eixo y.
# Usando o base plot system do R
plot(n_vector,MSE_vector,col=1, pch=1, main = 'Erro Quadrático Médio', 
     col.main = 'black', 
     xlab = 'Número de Amostras',
     ylab = 'MSE')
grid()

Usando o pacote ggplot2

# Criação de um dataframe a partir dos dados
data = data.frame(n_vector, MSE_vector) 

# Uso do pipe (tidyverse) e ggplot2
data %>% ggplot(aes(x = n_vector, y = MSE_vector)) + 
  geom_point() +
  xlab('Número de Amostras') + 
  ylab('MSE') + 
  ggtitle('Comportamento do MSE')

  1. Faça uma análise: o impacto com o aumento ou redução de n é significativo para o modelo? Qual a justificativa?

Sim é significativo, quanto menor o valor de n maior o MSE e quanto maior o valor de n menor o MSE.

Exercício Computacional VIII

Conceito explorado: Introdução à Análise de Séries Temporais

Esse exercício tem o objetivo de explorar o assunto de séries temporais de forma introdutória. Para isso, existem alguns pacotes que podem nos auxiliar no objetivo desse exercício: capturar séries temporais do mercado financeiro e realizar sua visualização. Abaixo, estão listados três pacotes relacionados ao mercado financeiro que podem ser instalados e carregados na linguagem R. • install.packages(“quantmod”) • install.packages(“xts”) • install.packages(“moments”)

Esses pacotes foram desenvolvidos exclusivamente para modelagem financeira quantitativa na lin- guagem R e permitem capturar séries temporais sobre as cotações de ações do mercado financeiro.

Especificamente, estude e utilize a função getSymbols do pacote “quantmod” para obter séries tem- porais de diversas empresas presentes na bolsa de valores a partir de uma janela de tempo fornecida.

Essa função consegue obter os dados diretamente das fontes “Yahoo Finance” (ainda ativo) e “Google Finance”, que disponibilizam os dados gratuitamente. Utilize o nome “yahoo” para designar a fonte de dados na função getSymbols.

## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## The following object is masked from 'package:leaflet':
## 
##     addLegend
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Version 0.4-0 included new data defaults. See ?getSymbols.
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
  1. Use a função getSymbols do pacote quantmod para capturar as cotações de ações da empresa Petrobras de janeiro/2020 até os dias atuais.
  dt_inicio = as.Date("2020-01-01")
  # Busca data atual do sistema.
  dt_fim    = Sys.Date() %>% as.Date()
  ?getSymbols
  
  getSymbols("PETR4.SA", src = "yahoo", from = dt_inicio, to = dt_fim, warnings = FALSE)
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
## Warning: PETR4.SA contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.
## [1] "PETR4.SA"
  # Verificando o tipo de objeto 
  class(PETR4.SA)
## [1] "xts" "zoo"
  is.xts(PETR4.SA)
## [1] TRUE

Mostrando o formato dos dados

summary(PETR4.SA)
##      Index            PETR4.SA.Open   PETR4.SA.High    PETR4.SA.Low  
##  Min.   :2020-01-02   Min.   :11.07   Min.   :12.18   Min.   :10.85  
##  1st Qu.:2020-02-07   1st Qu.:16.69   1st Qu.:17.09   1st Qu.:15.84  
##  Median :2020-03-18   Median :19.49   Median :19.66   Median :19.11  
##  Mean   :2020-03-18   Mean   :21.88   Mean   :22.29   Mean   :21.36  
##  3rd Qu.:2020-04-27   3rd Qu.:29.05   3rd Qu.:29.54   3rd Qu.:28.75  
##  Max.   :2020-06-03   Max.   :30.88   Max.   :31.24   Max.   :30.47  
##                       NA's   :1       NA's   :1       NA's   :1      
##  PETR4.SA.Close  PETR4.SA.Volume     PETR4.SA.Adjusted
##  Min.   :11.29   Min.   : 25397500   Min.   :11.29    
##  1st Qu.:16.44   1st Qu.: 48112975   1st Qu.:16.44    
##  Median :19.39   Median : 80764500   Median :19.39    
##  Mean   :21.81   Mean   : 88251612   Mean   :21.81    
##  3rd Qu.:29.15   3rd Qu.:116799625   3rd Qu.:29.15    
##  Max.   :30.81   Max.   :254813800   Max.   :30.81    
##  NA's   :1       NA's   :1           NA's   :1

Omitindo os NAS

petrobras = na.omit(PETR4.SA)
head(petrobras)
##            PETR4.SA.Open PETR4.SA.High PETR4.SA.Low PETR4.SA.Close
## 2020-01-02         30.51         30.70        30.31          30.70
## 2020-01-03         30.88         31.24        30.45          30.45
## 2020-01-06         30.43         30.94        29.95          30.81
## 2020-01-07         30.82         30.88        30.47          30.69
## 2020-01-08         30.69         30.77        30.24          30.50
## 2020-01-09         30.47         30.62        30.25          30.40
##            PETR4.SA.Volume PETR4.SA.Adjusted
## 2020-01-02        37774500          30.69834
## 2020-01-03        71595600          30.44835
## 2020-01-06        81844000          30.80833
## 2020-01-07        32822000          30.68834
## 2020-01-08        48215600          30.49835
## 2020-01-09        36102700          30.39835
# Observando o dataset - nós temos a cotação de abertura (open), a mais alta (high), mais baixa (low)
# a cotação de fechamento (close), finalizando com volume e cotação ajustada.
summary(petrobras)
##      Index            PETR4.SA.Open   PETR4.SA.High    PETR4.SA.Low  
##  Min.   :2020-01-02   Min.   :11.07   Min.   :12.18   Min.   :10.85  
##  1st Qu.:2020-02-06   1st Qu.:16.69   1st Qu.:17.09   1st Qu.:15.84  
##  Median :2020-03-18   Median :19.49   Median :19.66   Median :19.11  
##  Mean   :2020-03-18   Mean   :21.88   Mean   :22.29   Mean   :21.36  
##  3rd Qu.:2020-04-27   3rd Qu.:29.05   3rd Qu.:29.54   3rd Qu.:28.75  
##  Max.   :2020-06-03   Max.   :30.88   Max.   :31.24   Max.   :30.47  
##  PETR4.SA.Close  PETR4.SA.Volume     PETR4.SA.Adjusted
##  Min.   :11.29   Min.   : 25397500   Min.   :11.29    
##  1st Qu.:16.44   1st Qu.: 48112975   1st Qu.:16.44    
##  Median :19.39   Median : 80764500   Median :19.39    
##  Mean   :21.81   Mean   : 88251612   Mean   :21.81    
##  3rd Qu.:29.15   3rd Qu.:116799625   3rd Qu.:29.15    
##  Max.   :30.81   Max.   :254813800   Max.   :30.81
  1. Use a função candleChart do pacote quantmod para fazer a visualização da série temporal das cotações fechadas da Petrobras no período considerado. Pesquise e explore essa função; esse resultado é chamado de gráfico de velas (amplamente conhecido e usado em análises do mercado financeiro).
petrobras_fechamento <- petrobras$PETR4.SA.Close
?candleChart

candleChart(petrobras_fechamento, name = "Fechamento da ações PETR4.SA")

Gerando o candlechart pela função plot

plot(petrobras_fechamento, main = "Fechamento Diário - Ações da Petrobrás", 
     col = "red", xlab = "Tempo", ylab = "Preço da Ação", major.ticks = "months")

  1. Explore a série temporal obtida da empresa, teste outros períodos de tempo e identifique o significado dos campos trazidos da série.
# basic example of ohlc charts
df <- data.frame(Date=index(petrobras),coredata(petrobras))

fig <- df %>% plot_ly(x = ~Date, type="candlestick",
          open = ~PETR4.SA.Open, close = ~PETR4.SA.Close,
          high = ~PETR4.SA.High, low = ~PETR4.SA.Low) 
fig %>% layout(title = "Ações da Petrobrás")

Os dados contam com os valores de abertura/fechamento/máxima e mínima de ações da petrobrás.

  1. Use a função addBands, também do pacote quantmod, para plotar diretamente limites superior/inferior sobre a série temporal do item 1). Podemos parametrizar a função fornecendo:
  1. o período da média móvel e ii) a quantidade de desvio padrão relacionados com os limites.
# Candle com toda a série
candleChart(petrobras, name = "Evolução da ação PETR4.SA")

# Adição de limites superior/inferior sobre a série temporal
# Exemplo: bandas de bollinger
?addBBands
numero_desvios      = 1
periodo_media_movel = 20
addBBands(periodo_media_movel, numero_desvios)

https://www.bussoladoinvestidor.com.br/bandas-de-bollinger/

  1. A partir os itens anteriores, faça uma análise dos preços de cotações da Petrobras e verifique os reflexos do momento atual que estamos passando (março/abril/2020).

Pelo gráfico é possível observar um forte movimento de quedra no preço da ação, inclusive em alguns períodos ficando abaixo no nível inferior da banda de bollinger.